A computational model of stem cells’ decision-making mechanism to maintain tissue homeostasis and organization in the presence of stochasticity

The maintenance of multi-cellular developed tissue depends on the proper cell production rate to replace the cells destroyed by the programmed process of cell death. The stem cell is the main source of producing cells in a developed normal tissue. It makes the stem cell the lead role in the scene of a fully formed developed tissue to fulfill its proper functionality. By focusing on the impact of stochasticity, here, we propose a computational model to reveal the internal mechanism of a stem cell, which generates the right proportion of different types of specialized cells, distribute them into their right position, and in the presence of intercellular reactions, maintain the organized structure in a homeostatic state. The result demonstrates that the spatial pattern could be harassed by the population geometries. Besides, it clearly shows that our model with progenitor cells able to recover the stem cell presence could retrieve the initial pattern appropriately in the case of injury. One of the fascinating outcomes of this project is demonstrating the contradictory roles of stochasticity. It breaks the proper boundaries of the initial spatial pattern in the population. While, on the flip side of the coin, it is the exact factor that provides the demanded non-genetic diversity in the tissue. The remarkable characteristic of the introduced model as the stem cells’ internal mechanism is that it could control the overall behavior of the population without need for any external factors.

A computational model of stem cells' decision-making mechanism to maintain tissue homeostasis and organization in the presence of stochasticity

Najme Khorasani 1* & Mehdi Sadeghi 2
The maintenance of multi-cellular developed tissue depends on the proper cell production rate to replace the cells destroyed by the programmed process of cell death. The stem cell is the main source of producing cells in a developed normal tissue. It makes the stem cell the lead role in the scene of a fully formed developed tissue to fulfill its proper functionality. By focusing on the impact of stochasticity, here, we propose a computational model to reveal the internal mechanism of a stem cell, which generates the right proportion of different types of specialized cells, distribute them into their right position, and in the presence of intercellular reactions, maintain the organized structure in a homeostatic state. The result demonstrates that the spatial pattern could be harassed by the population geometries. Besides, it clearly shows that our model with progenitor cells able to recover the stem cell presence could retrieve the initial pattern appropriately in the case of injury. One of the fascinating outcomes of this project is demonstrating the contradictory roles of stochasticity. It breaks the proper boundaries of the initial spatial pattern in the population. While, on the flip side of the coin, it is the exact factor that provides the demanded non-genetic diversity in the tissue. The remarkable characteristic of the introduced model as the stem cells' internal mechanism is that it could control the overall behavior of the population without need for any external factors.
One of the major transitions in evolution is the major step from mere unicellularity into the new world of multicellularity. To understand this, one has to study two main characteristics of multicellular organisms that appeared during this transition, namely differentiation, and self-organization. These two processes are so tied up during development. To be more specific, through the self-organization process, genetically homogeneous differentiated groups of cells are organized into the tissues with different shapes and functionalities that can all together work as an overall patterned structure. For the proper functionality of the fully formed organism in a developed normal tissue, it is vital to have a strategy to generate the right proportion of different types of specialized cells, distribute them into their right position, and maintain the organized structure in a homeostatic state.
In a developed tissue, stem cells are the main source of producing cells 1 . Therefore the main focus of this project is understanding the stem cells' internal mechanism which controls their decisions. The decisions which satisfy the primary goal of life: keeping the living organisms alive. Stem cells are characterized by their capacity to self-renew and differentiate into more specialized cell types 2,3 . Any imbalance between the tissue demand and these two can lead to dysfunctional tissues or tumorigenesis. Hence, here we introduce a computational model to understand the underlying mechanisms that orchestrate the stem cell proliferation/differentiation balance to regulate the non-genetic diversity, as well as maintain the spatial pattern in the dynamic stochastic environment of a living tissue 4 . This will be a step toward the big picture of our final project: introducing a comprehensive model to understand the stem cells' internal activities regulating their final fate. That helps us study mysterious phenomena such as aging, and cancer, and manufacture artificial organs and develop stem-cell-based therapies in regenerative medicine.
That is for sure that, a non-genetic diversifying factor is required to produce hundreds of different cell types, in an organized pattern, from genetically similar stem cells. In this circumstance, stochasticity is advantageous and enormous number of specialized cells demanded in the tissue by going through many rounds of cell division 58 . Therefore, the production of thousands of differentiated cells could be maintained with a low rate of stem-cell division rate. This is beneficial in the sense that the larger number of stem cell division cycles, the greater is the risk of stem-cell mutations. This way, the presence of progenitor cells can prevent the accumulation of mutations which is sufficient to cause cancer. Therefore, it is assumed that the developing tissue that we are studying is made up of progenitor cells (P) as well as stem cells, and differentiated cell type 58,59 . In our model, the stem cells can self-renew and also can be divided into progenitor cells which present a bistable system. Progenitor cells can give rise to only one or a few types of specialized non-dividing cells. Without loss of generality, here the progenitor cells can divide into two non-dividing differentiated cells termed A, and B, which die after several days or weeks 59 . From the dynamical system point of view in our model, the progenitor cells are studied as a tristable system, which is also biologically reliable.
Investigations of tissue regeneration unveil a remarkable degree of flexibility, with progenitor cells able to recover the stem cell presence following injury, and even under normal conditions. It is also has been shown that reversible transfer of cells provides a more viable mechanism to maintain homeostasis under all sort of conditions [60][61][62][63][64][65][66][67][68][69][70] . That being the case, in this model stem and progenitor cells can switch stochastically between two states: stem cells could transit to progenitor cells and vice versa with a fixed rate. In addition, it is assumed that stem and progenitor cells diminish with rate of γ S , and γ P . Within the described framework, the model dynamics could be described as follows: where η , η s , and η p denote the rates of stem cell's three different division types, ω p,s represents the transition rates from S to P, and S to P, respectively, and p , A , B , µ d , µ A , and µ B denote the rates of progenitor cell's six different division types. The last two processes denote the rates, γ A , and γ B , at which A, and B cells commit to death.
The time evolution of the average densities of cell type S, cell type P, cell type A, and of cell type B-where n S , n P , n A , and n B are cell numbers normalized by volume-is given by In steady state: then, As n * S > 0 , n * P > 0 , and w S > 0 then: = n S η S − n S η P + n P w S − n S w P − n S γ S , ∂ t n P = g(n S , n P , n A , n B ) = n S η + 2n S η P + n S w P − n P w S − n P γ P − n P (− P + µ d + µ A + µ B ), ∂ t n A = h(n S , n P , n A , n B )  58 . In other words, after a limited number of cell divisions, they terminally differentiate to fulfill their responsibility in producing specialized cells in a normal mature tissue 59 . Consequently, it is interpreted that the progenitor cells' potential to self-renew lowered through generations. In our model, it is made to happen by setting the parameters in a way that the conditions in Eq. (11) are satisfied. Because, in contrast to the second condition in Eq. (10), when the conditions in Eq. (11) are satisfied, through generations it becomes more and more probable for a progenitor cell to be differentiated than self-renew. It can be easily interpreted as the limited capacity of proliferation in progenitor cells.
To study the stability of the model, it is needed to compute the Jacobian matrix of Eq. (6): where f (n S , n P , n A , n B ) , g(n S , n P , n A , n B ) , h(n S , n P , n A , n B ) , and q(n S , n P , n A , n B ) are denoted as F, G, H, and Q.
The system fixed point, (n * S , n * P , n * A , n * B ) , is stable if all four eigenvalues of J, or all the roots of equation det(J − I) = 0 , at the fixed point are negative 71 In this case, det(J − I) = (κξ − w s θ)υφ which is independent from the values of χ , and ψ . Therefore, one could say that from the stability point of view, our system in (6) is equivalent to two independent two-dimensional systems as follows: and In the former two-dimensional system, det(J ′ − I) = κξ − w s θ (first term in det(J − I)), and in the latter one, det(J ′′ − I) = υφ (second term in det(J − I)). Hence, if all four roots of these two equations are proven to be negative, the condition for the stability of the system in Eq. (6) is fulfilled.
Here, we aim to study the conditions which lead to the stability of the first system's fixed point, namely, (n * s , n * p ) . The time evolution of the average densities of cell type S, n s , cell type P, n P is given in Eq. (14). Clearly, the average total density of dividing cells (stem and progenitor cells), n, could be computed as follows: then, When we substitute n S = n − n P , the Eq. (14) becomes: where we introduced the functions (9) η S − η P − γ S − w P < 0.
   ∂ t n S = n S η S − n S η P + n P w S − n S w P − n S γ S , ∂ t n S = n S η S − n S η P + n P w S − n S w P − n S γ S , − n P (− P + µ d + µ A + µ B ), (18) ∂ t n P = f ′ (n P , n) ∂ t n = g ′ (n P , n), where, f ′ (n P , n) , and g ′ (n P , n) are denoted as F ′ , and G ′ . The fixed point (n * P , n * ) is stable if both eigenvalues of J ′ at the fixed point are negative. This is the case if For the system defined in Eqs. (18), and (19), the Jacobian reads as Eq. (25), where In steady state we could substitue w S = w * S (Eq. 8), and compute det(J ′ (n P , n)) , and tr(J ′ (n P , n)) as follows: It is clear that det(J ′ (n P , n)) > 0 , and tr(J ′ (n P , n)) < 0 as long as η ′ S (n) < 0 , and the conditions in Eqs. (9), and 11 are satisfied. Therefore, the fixed point (n * P , n * ) is stable. The Jacobian matrix of the second system (Eq. 15) reads: In this system, as det(J ′′ (n A , n B )) = γ A γ B > 0 , and tr(J ′ − (n A , n B )) = −γ A − γ B < 0 , the fixed point, ( n * A , n * B ) is stable. Therefore, it can be instantly concluded that the system fixed point of (n * S , n * P , n * A , n * B ) is stable.

Stem cells' internal mechanism.
The regulatory mechanism which provides the first part of the model dynamics in (1) is described by a set of ordinary differential equations (ODEs) which previously was used in several studies 7,32,40 . The following set of ODEs are employed to describe the stem cells' internal mechanism as a two-element bi-stable regulatory switch: In this model, It is assumed that the cell type is controlled by the relative amount of two cytoplasmic cell fate determinants, namely X s and Y s whose interactions can be described in a form of a bi-stable regulatory switch (see Fig. 1a). The dynamical behavior of the determinants X s and Y s is studied by considering their mutual repression effect which is modeled in the form of a Hill function 7,13,23 , and their degradation rate. Here, n is the Hill coefficient, β s is the effective rate of determinants synthesis, ι X s and ι Y s are inhibition rates, and γ is the degradation rate. The parameters of Eq. (26) are set in such a way that there will be two stable steady states, as shown in Fig. 1b, corresponding to two different cell fates, stem cell type S and progenitor cell types P. The number of determinants of X s ( Y s ) involved in attractor S (P) is much larger than those of Y s ( X s ). Figure 1b represents the domains of the two attractors, S, and P, with two different colors, yeloow and blue, respectively. Each daughter cell with a specific value of X s and Y s , right after birth, can be shown as a point in Fig. 1b. The value of X s and Y s determines to which attractor the cell will be absorbed, and based on that it defines the domains of two attractors. In other words, each cell fate can be determined and fixed exactly after division based on the number of determinants X s and Y s in the daughter cell.
Since, based on Eq. (26), ι x s , and ι y s are the parameters which determine the number of determinants of X s ( Y s ) involved in stem cells, it can be easily concluded that they are the parameters which control the rate of symmetric (19) f ′ (n, n P ) = n(η + 2η P + w P ) − → S + S , and S η P − → P + P ) directly. Furthermore, to have the first condition in (11) be satisfied in steady state, we set η S < η P . Therefore, for this purpose it is enough to set ι x s < ι y s . Besides, corresponding to η S = η s (n) , in our model ι x s is considered as a function of n, in a way that ι x s (n) ′ < 0.

Progenitor cells' internal mechanism.
The regulatory mechanism of the progenitor cells (shown in Fig. 1c) is described by the following set of ordinary differential equations (ODEs) as a two-element tristable system 7,32,40 : www.nature.com/scientificreports/ Likewise, It is assumed that the cell's final fate is determined by the relative amount of two cytoplasmic determinants, namely X p and Y p whose interactions can be described in a form of a tristable switch (see Fig. 1c). The dynamical behavior of the determinants X p and Y p is studied by considering their self-activation, and mutual repression effects which are modeled in the form of a Hill function 7,13,23 , and their degradation rate. In the same manner as the bi-stable system, n, β p , ι X s /Y s , and γ plays their role as the Hill coefficient, the effective rate of determinants synthesis, inhibition, and degradation rate, respectively. Besides, there are four more parameters, α x p , and α y p , as activation rates, and ε s 1 , and ε s 2 as "signalling effect coefficients" (will be discussed in this session). As shown in Fig. 1d, the parameters of Eq. (27) lead to a tristable steady-state system. In this system, there are three fixed points corresponding to three different cell fates, one progenitor cell type, P, and two differentiated cell types, namely, A, and B. The number of determinants X p , and Y p are in balance in P cells, where in A (B) cells the number of determinants X p ( Y p ) exceeds that of determinants Y p ( X p ). Figure 1d represents the domains of the three attractors, corresponding to P, A, and B cells, with three different colors, green, yellow, and blue respectively. As the number of determinants in the cell is updating, their corresponding trajectory in the phase plane is changing and finally reaches the domain of their attractor. It is worth mentioning that each cell fate can be finalized right after its birth based on the number of determinants X p and Y p in it.
Representing in Fig. 2a, it is assumed that there is a dish, occupied with four types of cells, S (cyan), P (green), A (yellow), and B (red), as the main scene for our system dynamics. As another assumption in our model, cell types A, and B produce signaling molecules namely, S 1 (yellow triangles), and S 2 (red triangles), respectively. Depicted in Eq. (27), theses signalling molecules could affect the self-activation rate of the progenitor cells through the parameters ε s 1 , and ε s 2 . The signaling effect coefficients read as follows: increases, ε s 1 ( ε s 2 ) increases in value. As a result, the self-activation effect of x p ( y p ) on itself grows, and the birth of A (B) daughter cells in that part of dish will be more probable. In the other words, by increasing the number of S 1 ( S 2 ) molecules, A (B) cells conquer the territory from other types of cell. Another point to mention is that a, and b values are set in such a way that territory of A (B) cells in dish remains unchanged.
As it is studied before 7 , β p is the parameter which determines the position of two boundaries in the force-field plane (in Fig. 1d). In other words, by shrinking (expanding) the progenitor cell's territory (the green domain), it controls the proportion of daughter cells which remain as progenitor cells, or differentiate to specialized cells. It can be easily concluded that it can directly control the rate of different types of divisions. To be specific, in our model, β p is set in such a way that the rate of symmetric division, P P − → P + P is much less than divisions in the forms of: P Therefore, the second condition in (11) will be satisfied in a steady state. To sum up, for each progenitor cell four reactions, X/Y synthesis/degradation, and one division process could occur, at each time step.
Signalling dynamics. As discussed in "Progenitor cells' internal mechanism", it is assumed that differentiated cells, A, and B secrete signaling molecules, S 1 , and S 2 respectively. The diffusion of signaling molecules between the pixels of the dish (Fig. 2a) can be governed by a set of reaction-diffusion equations as follows: Here D, and, k represent, respectively, diffusion coefficient and the rate of decay of the signaling molecules, and α s 1 ( α s 2 ) is the production rate. The number of signaling molecules at each point of the dish is shown by terms, s 1 , and s 2 , and their interactions are described in a form of a bi-stable regulatory switch (see Fig. 1e). The second term in Eq. (29) demonstrates the mutual inhibition effect of the signaling molecules, S 1 , and S 2 in the form of a Hill function. Here, n is the Hill coefficient, β is the effective rate of signaling molecules synthesis. Figure 1f represents the domains of the two attractors corresponding to S 1 , and S 2 , with two different colors, blue, and yellow, respectively. Each pixel of the dish with a specific value of S 1 , and S 2 can be shown as a point in Fig. 1f. The values of S 1 , and S 2 determines which attractor the cell will be absorbed to, and based on that it defines the domains of two attractors. Therefore, in the steady state, and in the deterministic environment, each pixel could only contain S 1 or S 2 , and not both.
In the simulations, the production of signaling molecules S 1 ( S 2 ) is proportional to the number of A (B) cells. When differentiated cells emerge in the dish, their corresponding signaling molecules diffuse in their environment where they interact with each other based on Eq. (29). The number of signaling molecules at any location in the dish determines how much a progenitor cell at that location is affected by the signal values. As shown in Eq. (27), the number S 1 ( S 2 ) will increase the birth rate of A (B) cells. Therefore, it is expected that if some parts of the dish are occupied by A (B) cells, it remains the same.
Gillespie algorithm. The time evolution of the system is captured by the Gillespie algorithm 7,13,44 which is known as the gold standard for simulating models whose stochasticity arises from the small discrete number of reactant molecules 72 . In each time step 24 different reactions can potentially happen. Table 1 demonstrates all the reactions and their corresponding propensity functions. In each iteration, one of the above-mentioned  www.nature.com/scientificreports/ processes occurs, time is updated, and the simulation continues to the point that all progenitor cells have gone through at least 50 divisions. The simulation starts with initializing the hypothetical dish (shown in Fig. 2a) with four types of cells and signalling molecules in an organized pattern. The number of determinants in stem (progenitor) cells, X S , and Y S ( X P , and Y P ), are initialized randomly from the corresponding attractor territory. The number of signalling molecules in each mesh is chosen corresponding to the initial pattern in the population. The meshes with A (B) cells contains maximum number, almost 100, of S 1 ( S 2 ) signalling molecules. In a mesh occupied with a stem (progenitor) cell, four reactions, production/degradation of determinant X S ( X P ), and production/degradation of determinant Y S ( Y P ) can potentially happen. In the deterministic manner, the ODE in Eq. (26) (Eq. 27) provides the exact description of these four reactions in our bistable (tristable) system. The propensity function of production/degradation reaction of determinant X S ( Y S ) is determined based on the positive/negative term in the first (second) ODE in Eq. (26) (rows 1-4 in Table 1). The propensity functions of production/degradation reaction of determinants X P , and Y P follow the same rule (see rows 9-12 captured based on Eq. (27)). At each time step, based on the probabilities corresponding to the propensity functions, the Gillespie algorithm determines which reaction occurs. As a result, in the current mesh, the determinants of the stem (progenitor) cell decrease or increase in number.
In a mesh with signalling molecules, based on Eq. (29), other than production/degradation of S 1 , and S 2 , diffusion is the fifth reaction to happen. The propensity functions for production/degredation of S 1 ( S 2 ), are determined based on second/third term in the first (second) ODE in Eq. (29) (see rows 19-22 in Table 1). The propensity function of the diffusion process is equal to D/h 2 (rows 23, and 24 in Table 1). The production/ degradation of signalling molecules occur in the system similar to the ones for dividing cells' determinants. However, when diffusion is the selected candidate to occur in the system, one of the neighboring meshes of the current mesh is selected, and the number of signalling molecules in the current mesh is decreased while the number of that in the neighboring mesh increases. The neighbor with much more number of signalling molecules comparing to the current mesh is more probable to be selected for the diffusion process. It is worth noting that, the above-mentioned propensity functions representing high order reactions could be used only as an approximation with Gillespie algorithm 73 .
The propensity functions corresponding to the rest of reactions are chosen as constant rates in such a way that satisfies the conditions introduced in "Mathematical model of the system". In the case of death, the cell in the current mesh is omitted from the population. In the case of division, one of the empty neighboring meshes is chosen for one of the newborn daughter cell. It is assumed that the distribution of determinants in each daughter cell is binomial 73 with parameters specified according to the whole number of determinants in the mother cell, and probability of success for each trial, p = 1/2 ( ∼ B(#X, 1/2) , or ∼ B(#Y , 1/2) , respectively).
All above-mentioned reactions are discussed before except the 8th, and 16th reactions, the stem, and progenitor cells' movement. In a real medium, when a cell divides its offspring could push their neighboring cells to make some space, and there is always a movement in the dish. However, reflecting these actions is beyond our simplified model. In our model we need this movement otherwise, in the regions where progenitors divide to non-dividing cells, the system would be blocked. In this case, A, and B cells can border progenitor cells, and since there will be no space for newborn cells, it restricts the cell division. To prevent a blocking system and to avoid physics complications to study the cell movements, it is simply assumed that dividing cells can move in the dish at a fixed rate. In the simplest form of the movement process, a dividing cell can change its position to a randomly chosen empty position in the dish.

Scoring algorithm.
To evaluate the strength of our model in maintaining the population pattern, we define two filters corresponding to the initial state of the medium, shown in Fig. 2b,c. These filters representing two matrices called template (corresponding to Fig. 2b) and penalty (corresponding to Fig. 2c) matrices from now on. The element values of the template matrix are equal to 1 in the middle disk of the current dish, they are equal to −1 in the outer ring, and the rest of the values are equal to 0. The elements of the penalty matrix are equal to 1 out of the current dish and are equal to 0 otherwise.
Each simulation starts with an initial state with a desired spatial pattern and continues to the point that all progenitor cells have gone through 50 divisions on average. We select 500 d × d shots out of all the states model meets during the simulation, and produce 500 2d × 2d matrices corresponding to them. Each pixel in the medium representation could be mapped to an element in a d × d matrix. Next, we add d/2 zeros to each of the four sides of the matrices. In all of these 500 matrices, elements corresponding to yellow/red pixels are set to 1/−1 , and the rest of the elements are set equal to 0. We slide the template, and penalty filters over each of the 500 matrices, multiply their corresponding values one by one and compute the summations. For each shot we do get two sets of d + 1 values, namely {t 1 , t 2 , . . . , t d+1 } , {p 1 , p 2 , . . . , p d+1 } , corresponding to template and penalty filters, respectively. The assigned score value for each shot is calculated as follows: Finally we normalize the score values in the range [0, 1].

Results
Maintaining the spatial pattern in the population. Aiming to study the capability of the model in maintaining the structural pattern of the population and in the presence of signalling molecules, the cells in the medium are initially organized in a circular pattern in which the dish is divided to two regions, a inner disk and an outer ring. It is assumed that the inner disk is mostly occupied with yellow A cells, where the outer ring is The first two simulations start with a population of n cells in a dish similar to the one presented in Fig. 2a. As shown in Fig. 3a, (and also Fig. 3d), the dish radius is initiated equal to 50 in the sense of placing at most 50 pixels on the dish radius. The inner disk radius is equal to 25. To compute the time evolution of the cell populations, a stochastic simulation, using Gillespie algorithm, is applied, and the simulation is terminated when the progenitor cells have gone through 50 divisions on average.
The first simulation is run in the absence of signaling molecules, and as expected the pattern could not be preserved in the population under this condition (Fig. 3b). However, in the second simulation, it is assumed that molecules A, and B produce signaling molecules as it was discussed before in "Signalling dynamics". The final state is shown in Fig. 3e. Visually studying the first (Fig. 3d) and last (Fig. 3e) states of the system, as well as the abundance of four cell types through time (shown in Fig. 3g), it is clear that in the second simulation, both abundance of cell types, and the initial pattern are maintained properly.
To evaluate the results by numbers, 500 shots are selected out of all the states model meets during the simulation. As explained in "Scoring algorithm", 500 corresponding scores are calculated and plotted in Fig. 3c,f, for the first and second simulations respectively. One could say that our scoring algorithm compares each of the 500 shots with the introduced filter in Fig. 2b. The filter which is chosen corresponding to the initial state of the system. In this case, diagrams in Fig. 3c,f demonstrates how much the population structural pattern differs from its initial pattern through the simulation. Figure 3c shows that at the end of the simulation the similarity between the initial and final states of our system is less than 10%. It clearly verifies that the initial pattern could not be maintained in the absence of intercellular interactions. Though Fig. 3f demonstrates the similarity score of ≃ 60% which indicates the model strength in maintaining the structural pattern in the presence of signaling molecules. Besides, in both cases, the score value diagram is saturated and it confirms the stability in the model which was promised in "Mathematical model of the system". To put it another way, our system reaches the point that the structural pattern and the number of different cells will not change anymore. Though the stochasticity is observable in all levels of the system dynamic.
Size of dish, and the inner disk affects the population spatial pattern. To show the key role of geometric confinement in maintaining spatial patterns in the population, we repeat the simulation with different sizes of inner disks as well as the whole dish. It is worth-mentioning that from now on all the simulations are done in the presence of signaling molecules. Figure 4a represents the results of five simulations with the same dish radius of 50, but different inner disk size, 5, 15, 25, 35, and 45 from left to right. The diagrams show that when the number of cells in one territory exceeds, they could easily invade the region occupied by the cells on the opposite side. To study the effect of the whole dish size in maintaining the organized pattern in the population, five simulations are run. Figure 4b demonstrates the results of these five simulations starting with five different dish radius, 10, 30, 50, 70, and 100, from left to right. In all of these simulations, the inner disk radius is www.nature.com/scientificreports/ always half of the dish radius. Figure 4b shows that there is a lower bound for the number of cells in the population, for cells to be able to defend their territory, and maintain the spatial pattern in the dish.
Model behaviour in the face of different population patterns. We simulated the model dynamics on other dish shapes and it is performed with any modification in order to study the effect of changing the geometry of the population. In Fig. 5a,b, the dish is initiated in the shape of a square, and rectangle, respectively. Studying the final results, one could say that there is an inward expansion of B cells at the corners of square and rectangular regions. Thus, the model behaves as it is expected based on experiments performed in previous researches 43 . Besides, in Fig. 5 it is clearly shown that by increasing the number of cells in the population, the inner region could be defended more easily by A cells. www.nature.com/scientificreports/ To challenge the model, we simulated the model in more complex patterns shown in Fig. 6. Here, again one fate has invaded the other fate in the corners. However, the model has been successfully capable of maintaining the initial pattern in the population.

Model in the face of injuries.
Studying the results, one could say that our model could properly maintain the organized pattern in the colony. To challenge the model, even more, a new experiment has been designed. In this experiment, in a state in which all progenitor cells have gone through 50 divisions on average(right-hand plot in the second row of Fig. 3), the cells in a part of the dish are diminished at once. This sudden cells death, Figure 5. The system behaviour for rectangle shape, and triangle shape dishes. The length of the middle area is always the half of the dish side length. The first, and second rows represent the initial, an final states of the system, respectively. The corresponding score diagram of five mentioned simulations are shown in the last row.  Fig. 7a,b). Studying the system behavior at this point, it is clear that the model could retrieve the initial pattern appropriately. However, one could say that the loss of cells in the boundaries of regions with different phenotypes is more challenging. Besides, as it is biologically expected the smaller size of the injury region, the more easily the system could recover.

Discussion
An organism's life depends on the precise functionality of its organs. On the other hand, the organs' proper functionality could not be achieved without the proper proportional tissue structure. Needless to say that there is a long way ahead to introduce a comprehensive model which controls the stem cell's decisions (stem cells as the main source of producing cells) to generate and also maintain the proper structure in the tissue. However, in this project, we took a step toward this tempting big picture and introduced a simple model which could properly maintain the pattern in the tissue.
Here, by focusing on the presence of stochasticity in all levels of cell activities, we define a regulatory switch as a mechanism to maintain both proliferation/differentiation balance and the spatial pattern in a hypothetical normal adult tissue. In the presence of progenitor cells an enormous number of specialized cells (to fulfill the tissue functionality) could be produced with a low rate of stem-cell division. It is pretty safe to conclude that it lowers the risk of mutation accumulation in stem cells! 58 . Besides, progenitor cells with the capability of reversible transfer between states destined for self-renewal or differentiation may provide a surprising degree of flexibility, with recovering stem cell population in the case of injury 60 . Therefore in the most simple model, it is assumed that this hypothetical developing tissue consists of progenitor cells other than stem cells and two differentiated cell types. As expected, stem cells could self-renew and differentiate to the intermediate progenitor cells, while progenitor cells could give birth to progenitor cells and two specialized cell types needed for the functionality of our hypothetical tissue.
The model is described in Eqs. (1), (2), (4), (5), in details. The computation in the rest of "Mathematical model of the system" proves that the described model could hit the homeostatic state under some reachable conditions. In other words, it is proven that starting with any initial condition, the desired proportion of specialized cells to satisfy the tissue functionality and stem, and progenitor cells as the main sources in the tissue, is reassuring. Besides, it shows that this condition could be maintained in a steady state, which is needed in any normal living tissue to survive 59 . Moreover, it reflects the flexibility of our model to describe any desired tissue with a different proportion of differentiated cells.
Two sets of ordinary differential equations are defined to describe the internal regulatory mechanism of the stem (Eq. 26), and progenitor cells (Eq. 27). As discussed before, although the regulatory networks could www.nature.com/scientificreports/ provide proliferation/differentiation balance and population heterogeneity and, as a result, reach and maintain a homeostasis state, they could not conserve the population structural organization. Therefore, the model is equipped to inter-cellular communication as the third issue required in order to maintain a spatial pattern. The results present evidence for the model's capability to maintain the organized spatial pattern in the population (Fig. 4a). The results also demonstrate that when one of the regions is much smaller than the other, the dominant cells could capture the territory of cells in the minority (the first and last plots in the second row of Fig. 4a,b). It is completely acceptable since as shown in Eq. (29), signaling molecules S 1 , and S 2 with symmetric mutual inhibition effects protect their territory. Therefore, It is necessary to break the symmetry in Eq. (29) to have a population containing regions with so much difference in population number. www.nature.com/scientificreports/ Here, studying different population sizes, consistent with results in previous studies 42, 43 the results suggest that dish size could influence maintaining spatial pattern in the population (Fig. 4b). Specifically, the smaller number of specialized cells in a dish, the fewer number of signaling molecules to defend their corresponding territory.
Importantly, we applied the model without any modification on other shapes (Figs. 5,6) to study the effect of changing the geometry of the colony. The fascinating outcome is that the model could easily maintain the spatial pattern in the population initiated in non-circular and even complex structures. However, an inward expansion of cells is observed at the sharp corners, which is in agreement with previous computational, and experimental predictions in previous studies 43 .
The death of cells caused by injury could cause organs failure due to an insufficient supply to fully restore tissue function. Therefore, as injuries are always a concern, shoulder to shoulder with the maintenance of the desired proportion of specialized cells and ordered spatial pattern, injury-induced repairability is demanded in the tissue 74 . Consequently, we computationally investigated the capability of our model to retrieve the initial pattern in the case of injuries. The results in Fig. 7 show that our model could successfully and without any further modifications restore the initial pattern, which is equivalent to restoring the tissue functionality.
The last but not least point is about the score diagrams shown in the last row of Figs. 4, 5, and 6. The scoring algorithm is designed to evaluate the results in numbers by representing how much the population spatial pattern differs from its initial pattern through 50 divisions. When the number of cells in two territories is both sufficient and in balance, the diagram decreases at the beginning, and at some point, it will get saturated. It clearly verifies that although the territories could not be kept in their exact initial pattern, the system reaches the point with the similarity score of 60% , and after that, the structural pattern and the number of different cells will not change anymore. It confirms the stability of the model. Here, a fundamental question could strike minds: what is the reason for having crooked boundaries instead of initial proper ones at the end of the simulations? The obvious answer is stochasticity. The stochasticity could be easily interpreted as a threat to the tissue functionality, which ties in with the tissue organized spatial pattern. On the other hand, stochasticity is the exact factor that brings demanded vital diversity out of the same genotypes in the tissue, the item, which facilitates the tissue functionality. These contradictory roles of stochasticity lead to one of the fascinating outcomes of this project: life is certainly indebted to a controlled amount of uncertainty.
In this model, by emphasizing the prominent role of stochasticity 5,7,13 in the developed adult tissue, we could maintain the spatial pattern to fulfill the tissue functionality other than the homeostasis state to have a cycling adult tissue 7,60 . We could even maintain more complex patterns and the results agreed with previous studies 43 . Besides, it has been mathematically proved that the introduced system could reach a steady state in all cases under some conditions. While injuries are always possible to happen 74 , with the presence of progenitor cells, we also demonstrate that our model could retrieve the right demanded proportion of different cells in their right position in the case of injuries. To summarize, we introduced a mechanism, by orchestrating the cell decisionmaking switch and in the presence of intercellular interactions could maintain the population's overall behavior with no need for any external factors.
Although it is vital to maintain the organized structural pattern in the tissue for the proper functionality, in a comprehensive model, it is also important to probe the stem cells' strategies to organize the specialized cells to a desired pattern in the colony. Thus, it could be interesting to investigate if our model could initiate pattern formation in the population and study the factors that could affect this process in future work. It could be a promising step toward the big picture of this project: generating organoids starting from stem cells.

Data availability
The software used to run all simulations was Python. The scripts and the data that support the findings of this study are available from the corresponding author upon request.